rm(list = ls())

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, fixest)

dat <- read_rds("data_not_upload/mz_1996_2016_clean.rds")

##

bw <- 10

dat <- dat %>%
  filter(year %in% c(2016)) %>%
  mutate(year = as.factor(year)) %>%
  filter(yob %in% seq(1971 - 9, 1971 + 9)) %>%
  mutate(group_var = case_when(
    east == 1 & male == 1 ~ "East German Men",
    east == 1 & male == 0 ~ "East German Women",
    east == 0 & male == 1 ~ "West German Men",
    TRUE ~ "West German Women"
  ))

# Create a new binary variable

dat$side_of_cutoff <- ifelse(dat$yob <= 1970.5, 0, 1)

# Plot

p1 <- ggplot(dat, aes(x = yob, y = abitur)) +
  stat_summary_bin(
    fun = mean, geom = "point", binwidth = 2,
    shape = 21, fill = "white"
  ) +
  geom_smooth(
    data = subset(dat, side_of_cutoff == 0), method = "lm", se = T, size = .5,
    col = "black"
  ) +
  geom_smooth(
    data = subset(dat, side_of_cutoff == 1), method = "lm", se = T, size = .5,
    col = "black"
  ) +
  facet_wrap(~group_var) +
  geom_vline(xintercept = 1970.5, linetype = "dotted") +
  theme_bw() +
  labs(
    x = "Birth year",
    y = "Share A-levels"
  ) +
  scale_x_continuous(breaks = seq(min(dat$yob), max(dat$yob), 2))

p1

## Parametric specification

dat$post <- ifelse(dat$yob >= 1971, 1, 0)
dat$yob_post <- dat$yob * dat$post

m1 <- feols(abitur ~ yob + post + yob_post,
  split = ~group_var,
  data = dat
)

etable(m1, keep = "yob_post")
